% Script written to extract the tide data from
% the tpxo7.2 model, multiply my Bz and write everithing to
% a file as an input to Alexei' code
% April 19, 2012

% The tides are provided as complex amplitudes of earth-relative
% sea-surface elevation for eight primary (M2, S2, N2, K2, K1, O1, P1, Q1),
% two long period (Mf,Mm) and 3 non-linear (M4, MS4, MN4) harmonic
% constituents, on a 1440x721, 1/4 degree resolution full global grid
clear

ncid = netcdf.open('/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_model_data/DATA/u_tpxo7.2.nc','NOWRITE');

% Variables
% 1 lon_u
% 2 lat_u
% 3 lon_v
% 4 lat_v
% 5 Ua
% 6 ua
% 7 up
% 8 Va
% 9 va
% 10 vp
% 11 URe
% 12 UIm
% 13 VRe
% 14 VIm

lon_u = netcdf.getVar(ncid, 1);
lat_u = netcdf.getVar(ncid, 2);
lon_v = netcdf.getVar(ncid, 3);
lat_v = netcdf.getVar(ncid, 4);
URe = netcdf.getVar(ncid, 11);
UIm = netcdf.getVar(ncid, 12);
VRe = netcdf.getVar(ncid, 13);
VIm = netcdf.getVar(ncid, 14);

% interploate the data to a common grid (the original grid is staggered)

latmat = lat_v;
lonmat = lon_u;

% remove the -90.125 values from the latmat

latmat(1,:) = [];
lonmat(1,:) = [];

%% Write flow data for each component to a file
% Dimention 1 is M2
% 1 M2,2 S2,3 N2,4 K2,5 K1,6 O1,7 P1,8 Q1
tidal_comp = ['M2';'S2';'N2';'K2';'K1';'O1';'P1';'Q1'];

for kk = 2:7,
    URe_inter = interp2(lon_u,lat_u, squeeze(URe(:,:,kk)), lonmat,latmat);
    UIm_inter = interp2(lon_u,lat_u, squeeze(UIm(:,:,kk)), lonmat,latmat);
    VRe_inter = interp2(lon_v,lat_v, squeeze(VRe(:,:,kk)), lonmat,latmat);
    VIm_inter = interp2(lon_v,lat_v, squeeze(VIm(:,:,kk)), lonmat,latmat);
    
    % Remove the NaNs at the edges
    
    VRe_inter(1:end,1) = VRe_inter(1:end,2);
    VIm_inter(1:end,1) = VIm_inter(1:end,2);
    
    % load the conductance map
    load /data/backup/mnair/conductivity/0.25x0.25_degree_resolution/0p25x0p25.cond;
    conduct = flipud(reshape(X0p25x0p25',[1440, 720])');
    
    % load the Bz
    % wmm evaluated at ellipsoid for 2007
    
    load /data/nobackup/mnair/wmmquarter/wmmquarter.txt
    
    Bz = (reshape(wmmquarter(:,5) , [1440,720]))';
    
    CUr = flipud(-VRe_inter.*Bz.*1e-9.*3.2); % Change sign since Br = -Bz
    CUi = flipud(VIm_inter.*Bz.*1e-9.*3.2); % Positive since tidal model has e^(i omega t) convention
    CVr = flipud(-URe_inter.*Bz.*1e-9.*3.2);
    CVi = flipud(UIm_inter.*Bz.*1e-9.*3.2);% Positive since tidal model has e^(i omega t) convention
    
    % test the data
    %imagesc_on_map(CUr,CUi,CVr,CVi);
    %CU = CU*1000; % Compression for e format 20 oct 2005
    %CV = CV*1000; % This multiplication is compensated by scale factor in source file
    %
    fid = fopen(['/data/nobackup/mnair/tides/tides_' sprintf('%s'...
        , tidal_comp(kk,:)) '_720_1440.source'],'wt');
    fid1 = fopen('/data/backup/mnair/ocean/header_deg0.25.txt','rt');
    
    %%the following format for alexei's x3dg is tested ok 20 oct 05
    for ii = 1:34,
        s=fgetl(fid1);
        fprintf(fid,'%s\n',s);
    end;
    s=fgetl(fid1);
    l = 1;
    for in = 1:720,
        for jn = 1:1440,
            fprintf(fid,'%12.4e',CUr(in,jn));
            if mod( l , 12) == 0,
                fprintf(fid,'\n');
            end;
            l = l + 1;
        end;
    end;
    fprintf(fid,'\n');
    for ii = 1:17,
        s=fgetl(fid1);
        fprintf(fid,'%s\n',s);
    end;
    s=fgetl(fid1);
    l = 1;
    for in = 1:720,
        for jn = 1:1440,
            fprintf(fid,'%12.4e',CUi(in,jn));
            if mod( l , 12) == 0,
                fprintf(fid,'\n');
            end;
            l = l + 1;
        end;
    end;
    fprintf(fid,'\n');
    for ii = 1:17,
        s=fgetl(fid1);
        fprintf(fid,'%s\n',s);
    end;
    s=fgetl(fid1);
    l = 1;
    for in = 1:720,
        for jn = 1:1440,
            fprintf(fid,'%12.4e',CVr(in,jn));
            if mod( l , 12) == 0,
                fprintf(fid,'\n');
            end;
            l = l + 1;
        end;
    end;
    fprintf(fid,'\n');
    for ii = 1:17,
        s=fgetl(fid1);
        fprintf(fid,'%s\n',s);
    end;
    s=fgetl(fid1);
    l = 1;
    for in = 1:720,
        for jn = 1:1440,
            fprintf(fid,'%12.4e',CVi(in,jn));
            if mod( l , 12) == 0,
                fprintf(fid,'\n');
            end;
            l = l + 1;
        end;
    end;
    fprintf(fid,'\n');
    for ii = 1:37,
        s=fgetl(fid1);
        fprintf(fid,'%s\n',s);
    end;
    s=fgetl(fid1);
    fclose(fid);
    fclose(fid1);
    
    fclose all;
    
end;

%% The codes below does the same thing as above, but at a spatial resolution 1x1 degree globally


% Script written to extract the tide data from
% the tpxo7.2 model, multiply my Bz and write everithing to
% a file as an input to Alexei' code
% April 19, 2012

% The tides are provided as complex amplitudes of earth-relative
% sea-surface elevation for eight primary (M2, S2, N2, K2, K1, O1, P1, Q1),
% two long period (Mf,Mm) and 3 non-linear (M4, MS4, MN4) harmonic
% constituents, on a 1440x721, 1/4 degree resolution full global grid

ncid = netcdf.open('/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_model_data/DATA/u_tpxo7.2.nc','NOWRITE');

% Variables
% 1 lon_u
% 2 lat_u
% 3 lon_v
% 4 lat_v
% 5 Ua
% 6 ua
% 7 up
% 8 Va
% 9 va
% 10 vp
% 11 URe
% 12 UIm
% 13 VRe
% 14 VIm

lon_u = netcdf.getVar(ncid, 1);
lat_u = netcdf.getVar(ncid, 2);
lon_v = netcdf.getVar(ncid, 3);
lat_v = netcdf.getVar(ncid, 4);
URe = netcdf.getVar(ncid, 11);
UIm = netcdf.getVar(ncid, 12);
VRe = netcdf.getVar(ncid, 13);
VIm = netcdf.getVar(ncid, 14);

% interploate the data to a common grid (the original grid is staggered)

latmat = repmat([-89.5:1:89.5]',[1,360]);
lonmat = repmat( [0.5:1:359.5],[180,1]);


%% Dimention 1 is M2
% 1 M2,2 S2,3 N2,4 K2,5 K1,6 O1,7 P1,8 Q1
tidal_comp = ['M2';'S2';'N2';'K2';'K1';'O1';'P1';'Q1'];

for kk = 1:7,
    URe_inter = interp2(lon_u,lat_u, squeeze(URe(:,:,kk)), lonmat,latmat);
    UIm_inter = interp2(lon_u,lat_u, squeeze(UIm(:,:,kk)), lonmat,latmat);
    VRe_inter = interp2(lon_v,lat_v, squeeze(VRe(:,:,kk)), lonmat,latmat);
    VIm_inter = interp2(lon_v,lat_v, squeeze(VIm(:,:,kk)), lonmat,latmat);
    
    % load the Bz
    % wmm evaluated at ellipsoid for 2007
    
    load /nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_model_data/wmmquarter/wmmquarter.txt;
    
    Bz = (reshape(wmmquarter(:,5) , [1440,720]))';
    
    bz_lat = repmat([-89.8750:0.25:89.875]',[1,1440]);
    bz_lon = repmat([0.125:0.25:359.875],[720,1]);
    % downsample Bz to 1x1 degree
    
    Bz_inter = interp2(bz_lon,bz_lat, Bz, lonmat,latmat);
    
    
    CUr = flipud(-VRe_inter.*Bz_inter.*1e-9.*3.2); % Change sign since Br = -Bz
    CUi = flipud(VIm_inter.*Bz_inter.*1e-9.*3.2); % Positive since tidal model has e^(i omega t) convention
    CVr = flipud(-URe_inter.*Bz_inter.*1e-9.*3.2);
    CVi = flipud(UIm_inter.*Bz_inter.*1e-9.*3.2);% Positive since tidal model has e^(i omega t) convention
    
    % test the data
    %imagesc_on_map(CUr,CUi,CVr,CVi);
    %CU = CU*1000; % Compression for e format 20 oct 2005
    %CV = CV*1000; % This multiplication is compensated by scale factor in source file
    %
    fid = fopen(['/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_model_data/sourcefiles/tides_' sprintf('%s'...
        , tidal_comp(kk,:)) '_180_360.source'],'wt');
    fid1 = fopen('/nfs/satmag_work/mnair/projects/ocean/OBEM/tpxo72_model_data/header.txt','rt');
    
    %%the following format for alexei's x3dg is tested ok 20 oct 05
    for ii = 1:34,
        s=fgetl(fid1);
        fprintf(fid,'%s\n',s);
    end;
    s=fgetl(fid1);
    l = 1;
    for in = 1:180,
        for jn = 1:360,
            fprintf(fid,'%12.4e',CUr(in,jn));
            if mod( l , 12) == 0,
                fprintf(fid,'\n');
            end;
            l = l + 1;
        end;
    end;
    fprintf(fid,'\n');
    for ii = 1:17,
        s=fgetl(fid1);
        fprintf(fid,'%s\n',s);
    end;
    s=fgetl(fid1);
    l = 1;
    for in = 1:180,
        for jn = 1:360,
            fprintf(fid,'%12.4e',CUi(in,jn));
            if mod( l , 12) == 0,
                fprintf(fid,'\n');
            end;
            l = l + 1;
        end;
    end;
    fprintf(fid,'\n');
    for ii = 1:17,
        s=fgetl(fid1);
        fprintf(fid,'%s\n',s);
    end;
    s=fgetl(fid1);
    l = 1;
    for in = 1:180,
        for jn = 1:360,
            fprintf(fid,'%12.4e',CVr(in,jn));
            if mod( l , 12) == 0,
                fprintf(fid,'\n');
            end;
            l = l + 1;
        end;
    end;
    fprintf(fid,'\n');
    for ii = 1:17,
        s=fgetl(fid1);
        fprintf(fid,'%s\n',s);
    end;
    s=fgetl(fid1);
    l = 1;
    for in = 1:180,
        for jn = 1:360,
            fprintf(fid,'%12.4e',CVi(in,jn));
            if mod( l , 12) == 0,
                fprintf(fid,'\n');
            end;
            l = l + 1;
        end;
    end;
    fprintf(fid,'\n');
    for ii = 1:37,
        s=fgetl(fid1);
        fprintf(fid,'%s\n',s);
    end;
    s=fgetl(fid1);
    fclose(fid);
    fclose(fid1);
    
    fclose all;
    
end;

%% Plotting contour maps of U and V as per Second reviwer's suggestion (Neesha paper)

